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-  two  decades,  much  research  effort  has  been 

.  devoted  to  the  application  of  the  stochastic  process  theory 
in  the  general  area  of  engineering  mechanics  and  structural 
engineering  for  the  purpose  of  predicting  the  dynamic  structural  \ 
performance  with  a  better  accuracy  and  of  assessing  the  over-  { 
all  structural  safety  with  a  better  reliability  by  considering 
more  realistic  analytical  models  of  load- structure  systems. 

Naming  only  a  few,  possible  applications  of  this  stochastic 
approach  include  (a)  analysis  of  panel  vibrations  of  aircraft 
and  submarines  induced  by  boundary  layer  turbulence,  (b)  analysis 
^  of  ship  oscillations  caused  by  ocean  waves,  particularly  during 
a  storm,  (c)  analysis  of  aircraft  response  to  gust  vertical 
velocity,  (d)  response  analysis  of  off-shore  structures  to  wave 
and  wind  forces,  (e)  statistical  strength  analysis  of  engineering 
materials,  particularly  of  modern  composite  materials,  with 
randomly  distributed  thermo-mechanical  properties,  (f)  analysis 
of  the  effect  of  randomness  in  geometrical  configuration  of  and 
mechanical  constraint  on  a  structural  component  due,  for  example, 
to  fabrication  errors  on  the  vibration  and  buckling  eigenvalues 
and,  (g)  study  of  random  surface  roughness  of  bridge  pavement 
-  and  airport  run-way  for  the  purposes  of  analyzing  the  vehicle  and 
aircraft  vibration  caused  by  the  roughness  and  of  stress  analysis 


of  pavement  systems  under  the  action  of  vehicles  and  aircraft^ 


^  ?  In  spite  of  the  recent  remarkable  advance  in  this  area  of  ^ 
study,  however,  the  present  state  of  art  still  leaves  a  number  ;  "  ; 
of  difficulties  unsolved  that  must  be  overcome  before  the 
approach  becomes  more  useful.  Such  problem  areas  include 
(!)  random  response  analysis  of  highly  nonlinear  structures, 

(2)  failure  analysis  of  structures  under  random  loading, 

(3)  analysis  of  extremely  complex  systems  and, 

(4)  random  eigenvalue  problems. 

The  recent  advent  of  high  speed  digital  computers,  however, 
has  made  it  not  only  possible  but  also  highly  practical  to  apply 
the  Monte  Carlo  techniques  to  a  large  variety  of  engineering 
problems.  The  present  paper  presents  a  technique  of  digital 
simulation  of  multivariate  and/or  multidimensional  Gaussian  random 
processes  (homogeneous  or  nonhomogeneous)  which  can  represent 
physical  processes  germane  to  structural  engineering.  The  paper 
also  describes  a  method  of  digital  simulation  of  envelope  functions 
Such  simulations  are  accomplished  in  terms  of  a  sum  of  cosine 
functions  with  random  phase  angles  and  used  as  the  basic  tool  in 
a  general  Monte  Carlo  method  of  solution  to  a  wide  class  of 
problems  in  structural  engineering,  particularly  those  mentioned 
above . 


2 


2.  Simulation  of  a  Random  Process:  A  Background 


A  basic  representation  of  a  homogeneous  Gaussian  (one¬ 


dimensional  and  one-variate)  random  process  fg  (x)  with  mean  zero  and 
spectral  density  (cu)  in  the  form  of  the  sum  of  the  cosine 
functions  has  existed  for  some  time  [1]; 


f  (x)  =  /2  2  cos  (cAj^x  -  §^) 
.  k=l 


where  §  are  random  angles  distributed  uniformly  between  0  and 
2rt  and  independent  of  (k  ^  j),  and 


Aj^  =  [S^  (cD^)  Am] 


=  (k  -  §)  Am 


with 


S^(m)  =  2So  (m)  (3) 

being  the  one-sided  spectral  density  function  (see  Fig.  1). 

For  digital  simulation  of  a  simple  function  f  (x)  of  f  (x)  and 
therefore  of  f^  (x) ,  Eq.  1  is  used  with  f  being  replaced  by 

Jv 

their  realized  values  cp  • 

iv 


f  (x)  =  /2  i:  Aj^  cos  (mij^x  -  cp^) 
k=l 


With  the  aid  of  Eqs.  1  and  4,  it  is  easy  to  show  that  Eq.  1  is 
ergodic  at  least  up  to  the  second  moment. 

Speaking  of  structure-related  applications,  however,  until 
Borgman  [2]  published  a  paper  simulating  ocean  surface  elevation 
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as  a  multidimensional  process  essentially  in  the  same  form,  little 
attention  had  been  paid  to  this  representation,  in  spite  of  its 
substantial  advantage  over  the  standard  method  in  which  a  random  ;  ■ 

process  was  digitally  generated  as  output  of  an  appropriate 
(analytical)  filter  subjected  to  a  simulated  white  process.  The 
use  of  this  filtering  technique,  although  limited  in  its  practical 
applications  to  a  one-variate  one-dimensional  process  has  dominated 
a  large  number  of  papers  involving  simulations  of  a  random  ^  > 

process. 

Practical  digital  simulation  of  a  multidimensional  process 
has  been  made  possible,  as  mentioned  above,  by  Borgman  [2]  in 
principle  through  the  preceding  expression  consisting  of  a  sum 
of  cosines  and  by  Shinozuka  [3]  whose  method  reduces,  in  case  of 
one-variate  one-dimensional  situations  to  the  use  of  the  following 
expression? 

-  2  4  ^ 

f(x)  =  aC:-)  i:  cos  (cD  X  -  cp  )  (5) 

^  k=l  ^  ^ 


where 


are  as  previously  defined  and 


are  realized  values 


of  random  frequencies  distributed  according  to  the  density 

2  . 

function  g  (cd)  =  (a3)/CT  with 


00 

=  J*  So  (oj)  dcD  (6) 

—  00 


Borgman  [2]  and  Shinozuka  [3]  also  investigated  the  digital 
simulation  of  the  multivariate  (but  one-dimensional)  process. 
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the  former  making  use  of  the  filtering  technique  and  the  latter  ' 
in  the  more  convenient  form  of  the  sum  of  the  cosine  functions.  - 

[4]  shows  that  the  autocorrelation  function 
R  (5;)  of  f  (x)  converges  to  Rj,  (f)  of  f^  (x)  In  the  form  of  1/N® 
as  N  -*  ®.  The  same  trend  in  convergence  can  be  shown  [4]  to  . 
exist  between  the  spectral  element  S  (cd)  of  f  (x)  and  Sq  (cd) 
of  fo(x)..' 

It  is  interesting  to  note  that  f  (x)  in  Eq.  1  can  be 
interpreted  as  a  canonical  expansion  of  a  Gaussian  process  f^  (x) 
with  mean  zero  and  spectral  density  (cu) .  To  see  this,  express 
fo  (x)  in  the  form  of  the  spectral  representation 


fo  (x)  =  J  e^“^  dZ  (cd) 


where  Z  (od)  ,  called  spectral  process,  is  orthogonal  in  the  sense 
that  the  increments  dZ  (cd^)  and  dZ(cD2)  are  uncorrelated  when 

“l  ^  ®2- 

Employing  the  orthogonal  condition  of  Z  (cd)  ,  the  auto¬ 
correlation  function  of  f.  (x)  is  found  to  be 


R,  (5)  =  E[fe  (x)  fg  (x+^)]  =  /  e^“^E|dZ(cD) 


where  E  indicates  the  expected  value. 

Assume  that  the  spectral  density  function  S,,  (cd)  exists. 

2 

Then  E|dZ(cD)  |  =  (cd)  dm,  and  Eq.  8  is  reduced  to  the  well- 

known  Wiener-Khintchine  relationship.  For  the  case  when  f„  (x) 
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is  :  real,-  Eq.  -7  becomes  ■ '[4] 

(x)  =  r  tcos  cuxdu  (oj)  +  sin  a)xdv(a»)]  (9) 

where  U  (co)  and  V  (cjo)  for  any  cu  s  0  are  two  mutually  orthogonal 
processes,  both  real  and  with  orthogonal  increments,  such  that 

=  E[dV(tD)]  =  S^Cco)  d(B 

^  out  that  if  one  defines 

dU(a3j^)=  [28^(03^)  Aco]^  cos  =  /2  cos 

"■  -(10)" 

dV  (cj^)  =  [28^  (oD^)  Acd]  ®  sin  =  /2  sin 

where  co  ,  Acc  and  I,  are  defined  in  Eqs.  1  and  2,  then  all  the 
k  .  k 

conditions  imposed  on  U  (co)  -  amd  V  (a>)  are  satisfied,  and  Eq.  1  ' 

is  basically  consistent  with  the  spectral  representation. 

It  is  seen  from  Eq.  9  that  a  homogeneous  process  is  addi- 
tively  built  up  by  orthogonal  oscillations  with  random  amplitudes. 

A  canonical  expression  of  a  real  random  process  f,  (x) 
can  be  written  in  the  form  [5] y 


fa  (x)  =  S  a  V  (x) 
k=i  ’= 


(11) 


where 

w^M 


a  are  uncorrelated  random  variables  with  mean  zero  and 
k 

are  real  (deterministic)  functions  of  x;  are  called 


the  coefficients  of  the  canonical  expansion  and  v,  (x)  the 


coordinate  functions. 
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multiplying  both  sides  of  Eq.  11  by  a^v  taking  the 
expected  values  and  using  the  fact  that  are  uncorrelated, 

one  can  show  that  .  -V.  V/ vv:./--*-- 


V  .  (x)  =  -  E[a  fo  (x)] 


where  D .  is  the  variance  of  a . ; 

1  .  3  . 


Dj  =  E[a^] 


At  this  point,  write  f^  (x)  in  Eq.  9  in  the  following 
approximate  form; 


f  (x)  =  E  [cos  05  X  du  (cD  )  +  sin  m  x  dv  (cu  )  ] 
-  -  K  K  JC  K 

k=l 


=  /2  E  (Aj^  cos  oj^x  cos  ^  si*'  si** 


/2  E  A^  cos  (o5^x  -  f^) 
k=l 


This  is  an  approximation  to  fo  (x)  since  the  integration  is 
replaced  by  a  summation  involving  a  finite  number  (N)  of  terms. 
The  degree  of  the  approximation  depends  on  (a)  whether  the 
cut-off  frequency  o)^  (see  Fig.  1)  is  large  enough  and  (b) 
whether  Aoa  in  Eq.  2  is  small  enough  so  that 


A.  =  f  S*  (od)  doj  =  S®  (<»  )  Ao) 

^  ‘(k-l)AcD  ^ 
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is  valid.  Being  the  sample  function  of  f  {x),  f  (xj  possesses  the 
same  approximate  nature  when  considered  as  the  sample  function 
of  :  f^  (x).  ■  ■  ■  ' 

within  the  context  of  this  approximation,  Eq.  14  can  be 
interpreted  as  an  (truncated)  canonical  expansion  of  £„  (x) 
in  the  form  of  Eq.  11  if  and  (x)  are  defined  as 

a  .  ,  =  /2  A,  cos  i  .  ,  ,  (x)  =  cos  od.x 

2d-1  :  3  .  23-1  3 

(j=l,2,...,N)  (16) 


a  =/2A.  sin$.  ,  V  .  (x)  =sincD.x 

2d  D  3  2d  D 


In  fact,  the  coordinate  functions  in  Eq.  16  follow  from 

Eq.  12  if  the  coefficients  a  in  Eq.  16  and  f  (x)  in  the  form 

Ic 

of  Eq.  14  are  used  therein.' 

It  can  then  be  shown  [5]  (again  within  the  context  of  the 
approximation  mentioned  above)  that  if  these  defined  in 

Eq.  16  are  used,  then  for  any  given  number  of  N  and  for  the 
particular  selection  of  a  in  Eq.  16,  Eq.  14  gives  the  best 

approximation  to  the  random  function  f,,  (x)  in  the  sense  that, 
among  all  possible  series  expansion  of  f^  (x)  having  the  same 
number  of  terms,  Eq.  14  minimizes  the  expectation  of  the  square 
of  the  residual  term  at  any  value  of  x. 

There  is  another  constraint  to  be  imposed  upon  Aco.  This 

constraint  is  due  to  the  periodicity  of  the  sample  function  f  (x) . 
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Obviously,  the  period  of  f  (x)  is  To  =  2rt/A<0  ;  and  therefore,^^^^^^'^^-^^^^^ 
depending  on  the  purpose  of  the-  simulation.  Am  has  to  be  so 
chosen  that  To  is  long  enough  for  that  purpose.  '  I  V-.> - 

A  significant  improvement  in  the  efficiency  of  digital 

simulation  has  recently  been  suggested  by  Yang  [6]  writing 

ReF(x).  (17)' ,  : 


in  which  Re P  (x)  represents  the  real  part  of  F (x)  and 


(x)  =1:  {[2  S“  (cc  )] 
k=l' 


-icp  im  t 
■  k-i  k 
!  ^e 


JC  ’ 

is  the  finite  complex  Fourier  transform  of  [2S^(mj^)]^  e 
with  and  cp^  defined  in  Eqs.  1  and  2.  The  advantage  of 
Eq.  18  is  such  that  the  function  F  (x)  can  readily  be  computed 
by  applying  the  fast  Fourier  transform  (FFT)  algorithm,  hence 
avoiding  the  time-consuming  computation  of  a  large  number  of 
cosine  functions. 

In  the  preceding  discussion,  the  spacing  Am  in  the  frequency 
domain  has  been  taken  as  constant.  This,  however,  does  not 
necessarily  have  to  be  always  observed.  It  is  possible  and  in 
fact  may  even  be  advisable  to  use  variable  spacings  depending  on 
the  extent  of  fluctuation  of  the  spectral  density  to  optimize 
the  number  of  cosine  terms  in  the  summation  in  Eq.  1;  finer 
spacings  in  those  domains  of  frequency  where  the  fluctuation  of 
the  spectral  density  is  more  rapid  and  coarser  spacings  elsewhere. 
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It  is  also  likely  that  such  variable  spacings  will-  increase  s  the 
period  of  the  sitnulated  process,  if  this'  is  done,  however 

the  simulation  through  the  FPT  technique  does  not  appear  to 
be  possible.  ' v/ ■  ■ 

same  reference  [6]  ,  Yang  also  proposed  to  simulate 
the  envelope  process  (t)  of  a  random  process  fo  (x)  by 
following  the  definition  of  the  envelope  process  [7] 

(t)  =  [fj  (x)  +  (x) 

where  (x)  is  the  Hilbert  transform  of  f^  (x)  and,  in  the 
present  case,  can  be  written  as 


fo  (x)  =  J  [sintut  dv (cd)  -  cos  cut  dv (cu) 3  (20) 

O 

/V  A 

It  then  follows  that  f^  (x)  can  be  simulated  as  f (x) y 

..  -N  ■  '■ 


.  i  (x)  = 

,/2  Z  A, 
k=l  ^ 

sin  (<u  X  -  cp,  ) 
k  k 

(21) 

Hence,  the  envelope 

process 

Vo  (t)  can  be  simulated  as 

V(t)  y 

2 

^2  4 

(22) 

v(t) 

=  [f  (x) 

+  f  (x)]^ 

As  an  example,  consider  the  response  process  y^  (t)  of 
single  degree  of  freedom  system  to  a  Gaussian  white  noise 
excitation  x^  (t)  with  the  constant  spectral  density  ; 

Yo  (t)  +  2(;cu„  y„  (t)  +  cuf  y,  (t)  =  x^  (t)  (23) 


It  is  well  known  that  the  spectral  density  of  y^  (t)  is 


-  2  2  2  2  2  2 
(co  “  CUo  )  +  AQ  CUo<U 
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and  the  standard  deviation  of  (t)  is 


JfSp  N.i 


(25) 


With  the  aid  of  Eq?-1,  21  and  22,  a  segment  of  sample 
function  y  (t)  of  simulated  process  y  (t)  for  y^,  (t)  and 
that  of  simulation  V(t)  of  the  envelope  process  (t)  are 
computed  and  shown  in  Pig.  2. and  3.  Fig.  2  is  for  the  case 

vhere  the  damping  coefficient  Q  =0.02  and  hence  the  process 
Yo  (t)  is  narrow-bando  This  fact  is  well  demonstrated  by  the 
smooth  behavior  of  the  sample  envelope.  When  the  local 
maxima  of  the  sample  envelope  do  not  coincide  with  the  local 
maxima  (peaks)  of  the  process,  they  reflect  the  local  minima 

(troughs)  of  the  simulated  process.  Fig,  3  shows  the  sample  functions  of 
V (t)  and  y(t)  for  q  =  0,5,  In  this  case,  the  process 
y^  (t)  is  substantially  wide-band  and  this  fact  is  clearly  seen 
from  the  much  wilder  fluctuation  of  both  simulated  envelope  and 
simulated  process,  although  the  simulated  envelope  surprisingly 
well  reflects  peaks  and  troughs  of  the  simulated  process  even  thou^ 
the  process  y^  (t)  is  wide-band.  In  terms  of  simulation 
efficiency,  the  computer  time  will  be  significantly  reduced  if 
one  is  interested  in  peak-  and  trough-values  of  the  process  and 
if  the  process  is  narrow-band,  since  then  the  smooth  nature  of 
the  envelope  function  makes  it  possible  to  use  much  larger 
interval  between  successive  time  instants  at  which  the  values  of 


th6  simulated  process  is  evaluated.  The  order  of  magnitude  of 
this  interval  can  be  that  of  the  apparent  period  of  the  process, 
which  obviously  is  much  too  large  for  simulation  of  the  process  : 
itself.  ' 

-  '  '  . ’-V.  .''-..''.-i-v.- 

^  write  f  (t)  as 

Vv\/V- ■ f  ImF  (t)  ■■ 

and  hence  Vi(t)  can  also  be  simulated  through  the  FFT  technique. 

It  appears  at  this  time  that  the  method  of  simulation 
considered  herein  (Eq.  1)  has  a  difficulty  in  achieving  a 
reliable  evaluation  of  the  first  passage  time  distribution 
when  the  threshold  value  is  much  larger  than  the  standard  de¬ 
viation  of  the  process.  Lyon's  work  [8]  points  to  this  fact, 
although  this  difficulty  is  by  no  means  unique  to  the  proposed 
method.  It  is  suggested  however,  that  a  further  investigation 
be  performed  on  this  point. 

The  Gaussian  property  of  the  simulated  process  (Eq.  1) 
comes  from  the  central  limit  theorem  because  it  consists  of  a 
sum  of  a  large  number  of  independent  functions  of  time  (see 
pp.  182  -  183  in  [1]).  Efficient  simulation,  or  straight  forward 
simulation  if  not  efficient^  of  a  non-Gaussian  process  appears 
to  be  an  open  problem  at  this  time  unless  the  process  is  restricted 
to  a  certain  class  of  processes  such  as  the  filtered  Poisson 
process. 

In  the  following  sections,  a  method  for  digital  simulation 
of  multidimensional  and/or  multivariate  processes  are  briefly 
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described.  However,  for  these  cases,  the  rigorous  discussions 
on  the  interpretation  as  canonical  expansion,  the  use  of  the 
PFT  technique  in  actual  digital  computation,  the  envelope 
simulation,  the  problem  of  the  first  excursion  time,  the  simu-  ‘  ^ 

lation  of  non-Gaussian  processes,  the  convergence  of  the  auto¬ 
correlation  function  and  the  spectral  density  of  the  simulated 
process  to  the  respective  target  values,  etc.,  are  mostly  subjects 
of  future  studies. 
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3.  Simulation  of  a  Multidimensional  Hoinoaeneous  Process 

autocorrelation  function  of  an  n^ditnensional  homogeneous 
real  process  f^  (x)  defined  by 

.;,Ro'  ^)V=  E[fo'.t?i),  ■■'■■  -V 

is  even  in  f  (symmetric  with  respect  to  the  origin  of  the 
n-dimensional  space) 

(z)  =  R„  (-^)‘  (27) 

where  x.  and  x_  are  space  vectors  and  C  =  “  x  is  the 

-1  -2  ^  2  -2  -1 

separation  vector.  Assume  that  the  n- fold  Fourier  transform 
of  Rq  (p)  exists.  The  spectral  density  function  of  f,,  (x)  is 
then  defined  as 

S,  (ffi)  = - /  R,  (5)  -dg  (28) 

■  (2n)"  -  .  .  ~ 

where  _co  is  the  frequency  (wave  number)  vector  and  ^  •  e;  is 
the  inner  product  of  q  and  z  ,  and,  for  simplicity 

/  (  )  <3l  =  J*  .  /  (  )  df^d^2-“‘^5n 

mmCD  *■■■00  <•■■00 

with  n  being  the  dimension  of  the  vector  It  follows  from 

Eq.  27  that 

00 

/  R„  (^)  sin(cD  •  ^)  =  0 

—  00 

and,  therefore,  from  Eq.  28 

S„  (cb)  =  So  (-cd)  (29) 
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Then 


■■  '\t:  >>^v 


So  (iS)  = 


(2  k) 


n 


/  Ro  {?:)  cos(a5  •  5)  df  ■  '  (3  0) 


and  is  real. 


can  be  shown  (9]  that  (§)  is  nonnegative  definite 
and  therefore  it  has  a  nonnegative  n-fold  Fourier  transform; 


So  (m)  5  0 


(31) 


Based  on  these  properties  of  Sq  (o?)  /  a  method  of  simulating 
fo  (x)  is  proposed  in  the  following; 

Consider  an  n-dimensional  homogeneous  process  with  mean 
zero  and  spectral  density  function  Sg(^)  which  is  of  insigni¬ 
ficant  magnitude  outside  the  region  defined  by 

<  oi.  i*  .a:  ^  05  <  <*> 

—  — u. 

and  denote  the  interval  vector  by 

.0)  “"03  .  CD  —CD  CD  ..~CD 

(Acd^,  Aa5^)  =  (  . . ^  ^  nf  ^  ^32^ 

1  2  n 

where  usually  05  =  ,  Then  the  process  can  be  simulated  by 

i/  u 


the  series 


N-  N  N 
1  2  n 


f  (x)  =  /2  E  E  ...  E  [So  (05  -u  ^  Acd.  Ao)  .  .  .  Am  ] 

-LK_  ZK^  nK  J.  2  n 


^1  ^  ^^2  ^  ^n~^ 


1  2 


n 


cos  ...  ^  ,  )  (33) 

1  2  n  1  2  n 
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where 


k  k  ~  independent  random  phase  uniformly 
:  2*  *  * ,  n, 

distributed  between  0  and  2^ 


“ik.  “  “ii'"  ''“i 


k^  =  1, 2 , . .  .  ,N^  i  =  1, 2, , . .  ,n 


As  in  the  one-dimensional  case,  the  digital  simu¬ 
lation  f  (x)  of  f  (x)  can  be  achieved  by  using  Eq.  33  with 

k-k  .  ..k  being  replaced  by  their  realized  values  cp  ,  ,  , 

1  “  n  k_  k  • • • k 

1  2  n 

To  avoid  the  lengthy  expressions  in  the  subsequent  discussion, 
f  (x)  will  be  written  in  the  following  compact  form: 


where 


f  (x)  =  /2  E  A(cd  )  cos  (ffi  •  X  +  cp.  ) 
k=l  ^  ^  ^ 


N  =  N,N  ...N 
12  n 


A(ai|^)  =  [So  (cD^)  =  [S,  (ci^)  ^  -  (35) 

It  is  noted  that  if  the  symmetric  condition  of  S,,  (^)  is  used, 

N  in  Eq.  33  can  be  reduced  by  one-half.  Furthermore,  if  the 

process  is  isotropic,  N  is  reduced  to  —  .  Fig.  4  illustrates 

2^ 

the  significance  of  A(m  )  for  two-dimensional  cases  where, 

JC 

however.  A,  ,  is  written  for  A(cd  m  ). 

It  can  be  shown  [4]  that  the  ensemble  average  of  f  (x)  is 
zero,  and  the  autocorrelation  fvinction  R(5)  of  f  (x) ,  becomes 


R(J^)  =  ^E^  A  (^)  cos(a5j^  .  §) 
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Upon  substituting  (^j^)  =  Sq  and  taking  the  limit  as 

N  -♦  «  (in  the  sense  that  N,  ,N- , . .  .  *  simultaneously);  one 

obtains  •  .  .  , -  -■ 


J  So  (cd)  cos  (m  •  _|)  dm  =  Ro  (_|) 


where  it  is  assumed  S„  (m)  =  0  for  m  <  cd.  and  m  >  m  . 

.  '~Jo  ^~u 

This  indicates  that,  when  the  ensemble  average  is  considered 
the  simulated  process  f(x)  possesses  the  target  autocorrelation 
Ro  (§)  and  therefore  the  target  spectral  density  S,  (m) . 

It  can  also  be  shown  [4]  that  the  temporal  (or  spatial)  mean 

<  f  (x)  >  is  zero  and  the  temporal  autocorrelation 
R*  (|)  =  <  f  (x)  f  (x  +  5)  >  becomes 


R*  (§)  =  E  A  (ax  )  cos  (m 
k=l  . 


As  N  -*  Eq.  38  becomes 

®  ' 

R*  (5;)  =  /  So  (m)  cos  (m  •  £)  dm  =  R,,  (^)  (39) 

**  —00 

From  Eqs.  36  and  38,  it  is  seen  that  the  process  f  (x)  in 
Eq.  33  is  ergodic  regardless  of  the  size  of  N.  This  makes  the 
method  directly  applicable  to  a  time  domain  analysis  in  which  the 
ensemble  average  can  be  evaluated  in  terms  of  the  temporal 
average. 

Note  that  the  simulated  process  is  Gaussian  by  virtue  of  the 
central  limit  theorem. 
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where  t  and  x  represent  the  time  and  the  distance  respectively 
and,  correspondingly,  a>  and  k  the  frequency  and  the  wave 
number.  It  is  known  that  such  process  (t,x)  is  a  satis¬ 

factory  model  of  a  fluctuating  part  of  wind  velocity  along  a 


straight  line  direction  x.  In  the  wind  study,  however,  it  is 
customary  to  consider  the  Fourier  transform  (cd,  ?)  of  the 
autocorrelation  (t,^)  =  E[fo  {t,x)  f^  (t+T,  x+  ^)]  only  with 
respect  to  t.  In  fact,  the  .form  of  (oa,  5)  consistent  with 


Eq.  40  is 


So  (03,  §)  -  2  *  2  2  4/3 

2ir  (1+  c  03  ) 


-a  I  031  I  ^ 


a  familiar  form  for  a  fluctuating  part  of  wind  velocity  at  the 
.  reference  altitude  of  33  feet  where  L  =  4000  ft.,  K  =  surface 
drag  coefficient,  a  =  constant,  C  =  1,/  with  being 

the  mean  wind  velocity  at  the  reference  altitude.  For  =40  mph^ 

a  =  0.02  ft*  sec  and  K  =  0.03,  the  sample  functions  f  (t,|)  of 

fo  (t, 0  are  computed  and  shown  in  Fig.  5  at  §  =0,  50  and  200  ft. 
One  can  easily  see  in  this  example  that  the  correlation  almost 
disappears  as  the  separation  §  increases  to  200  ft. 
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Consider  a  set  of  ra  homogeneous  Gaussian  n-dimensional  pro¬ 
cesses  fj(x)  ( j  =  1,2 , ....  . .  m)  with  meaih  zero  and  with  the  cross- 
spectral  density  matrix  S  (oi)  defined  by 


Sot  (ilj)  (co)  •••  (cn) 


Nina's’  ••• 


Where  s!,  (to)  is  the  n-fold  Fourier  transform  of  the  cross 
■  Dk  -  . 

correlation 

Due  to  the  fact  that  R*^(^)  =  one  obtains 


®jk<a>  -  \V<a) 


(43) 


where  the  bar  indicates  the  complex  conjugate.  __ 

The  matrix  s“(m)  is  therefore  Hermitian.  As  in  the  case  of 
a  one-dimensional  multivariate  process [10} , it  can  be  shown  [4] 
that  the  matrix  s“(<n)  is  also  nonnegative  definite. 

Suppose  one  can  find  a  matrix  H(co)  which  possesses 
n-dimensional  Fourier  transform  and  satisfies  the  equation 

S*(m)  =  H(m)  (44) 

where  S®  (m)  is  the  specified  target  cross-spectral  matrix  and  T 
indicate  the  transpose.  Then  f .  (x)  (j=l,  2,... ,m)  can  be 
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simulated  by  the  following  filtering  technique  [2,  li]  ; 

"i  V:  (45) 


m 


f.cx)  =  s 

3  lc  =  l  '’-oa  "■ 


where  h (x)  is  the  n-dimensional  Fourier  tj^onsofrm  of  ; 


■  ■  ■  ‘  **00 

and  'nj^(2.)  is  an  independent  n-dimensional  normalized  white 
noise  component  such  that 

E[Tlj  (Xj^)'nj^(x2)]  =  6(Xj^  “  2.2^^j]c 

With.'' 

6(x^‘-X2^  =  ^  ^^12”^22^'**'**®  ^^In  “  ^2n^ 

It  can  be  easily  verified  that  the  f ^ (x)  (j  =  l/...^m), 

as  simulated  by  Eq.  45,  satisfy  Eq.  44  and  thus  represent  the 
target  processes. 


To  find 

the 

matrix  H (m) 

in  an 

efficient 

way 

assume  that 

H(co) 

is 

a  lower 

triangular  matrix; 

(cd) 

0 

0  .  . 

.  0 

H(ai)  = 

•  •  • 

.  0 

\2  '®> 

•  •  • 

H 

ram 

Substituting  above  into  Eq.  44,  solutions  are  obtained  (see 
Ref , 12  for  similar  derivation)  as 
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is  tlie  determinant  of  a  submatrix  obtained  by  deleting  all  elements 
except  the  (1, 2, . .  .k-1,  j)-th  row  and  (1, 2,... k-l,]c)-th  column 
of  (cd)  . 

It  is  noted  that  the  above  decomposition  is  valid  only 
when  the  matrix  S°  (oi)  is  hermitian  and  positive  definite  as 
can  be  seen  from  Eq.  46. 

Because  the  cross-spectral  density  matrix  S° (m)  is 

known  to  be  only  nonnegative  definite,  special  consideration 

o 

is  needed  in  those  cases  where  S (^)  has  a  zero  principal 
minor.  For  the  discussion  on  this  point,  the  reader  is 
referred  to  Ref.  [4]. 
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Since  tJie  real  and  the  imaginary  parts  of  cross-spectral 
density  functions  are  respectively  even  and  odd  in  o),  it  can  be 
sbown  from  successive  substitutions  using  Eqs.  46  3nd47' 


Re  H.,  (to)  =  Re  H.,  (-to) 
3k  “  3k  ~ 

Im  H..  (to)  =-Ira  H.,  (-to) 
^jk^  3k  - 


for  j  >  k,  and 

H,  .  (to)  =  H.  .  (-to)  s  0 
33  -  33  -  • 

from  which  it  follows  that  h .  (x)  is  real 

3^ 


If  v/ritten  in  polar  form; 


then,  due  to  Eq.  48,  the  argument  ®  (s)  anti-symmetric 

in  m  ,  that  is 

with  © .  .  (to)  =  0  . 

33 

Once  H  (m)  is  computed  using  Eqs.  46  and  47,  then  instead 
of  passing  a  white  noise  vector  through  filters,  the  process 


f  j  (2?)  can  be  simulated  in  terms  of  the  following  series 
3  N 

(x)  =  E  E  |h.  cos[m  •  x  +  0.^(^.)+  f_J  (51) 


»  im 

1=1  i=i  ^  ^ 


"tag 


where  Am,  N,  and  are  essentially  the  same  as  defined 

previously  for  n-dimensional  processes  and 
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"iwH ;  (ci  j 


Q  ,  (m  )  —  tan  (  — - ^n'\  j 

3m'— i  \  Re  H.  i^l)  y  :  w-.'-y:  K',:y:y.^yfiy/. 


It  can  be  shown  [4]  that  the  processes  f ^  (x)  (j  -  l,  ****!^)* 

as  simulated  by  Eq.  51,  possess  the  target  cross-correlation 

functions  and  hence  the  target  cross-spectral  density,  with 

respect  to  an  ensemble  average. 

For  digital  simulation  of  sample  functions  of  f ^  (x) , 

Eq.  51  is  used  with  $  .  being  replaced  by  their  realized  values 

niK 
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simulation  of  Multidimensional  Nonhotnogeneous  Processes 


Simulations  of  nonstationary  processes  have  been  studied 
dealing  mostly  with  earthquake  ground  motion.  The  common  feature 
of  these  studies  is  that  a  nonstationary  process  is  simulated 
by  multiplying  by  an  envelope  function  the  stationary  process 
generated  either  by  filtering  a  white  noise  [13,14]  or  by  a 
series  of  oscillations  with  random  frequency  and  random  phase 


■[15,16,171. 

The  efficient  method  of  simulation  that  has  been  proposed 
for  multidimensional  homogeneous  processes  can  be  directly 
generalized  to  a  nonhomogeneous  process  characterized  by  an 
evolutionary  power  spectrum  as  introduced  by  Priestley  [18,19], 

It  was  seen  from  Eq.  9  that  a  homogeneous  process  is  addi- 
tively  built  up  by  orthogonal  oscillations  with  random  amplitudes. 
This  concept  of  orthogonal  components  can  be  extended  to  that  of 
the  evolutionary  process  f°  (x)  expressed  as 


CO 

f°  (x)  =  r  B(x,cu)[cos  cjDxdU  (cd)  +  sin  cDxdV(cu)]  (53) 

e  "  0 


where  B(x,a>)  is  a  deterministic  modulating  function  characterizing 
the  "nonhomogeneity "  of  the  process,  and  U  (a>)  and  V  (cd)  are  the 
same  as  defined  in  Eq.  10. 

Using  the  orthogonal  conditions  of  U  (cd)  and  V  (co) ,  the  mean 

square  of  f®  (x)  is  found  to  be 
© 


00 


E[f“(x)].^=  J*  (x,©)  S°  (cu)  do)  =  ,7  s°  (x,’^)  ,dw' 

,  where  S^(x,ai)  =  B  (x,cd)  8^(03)  is  defined  as  the  evolutionary,  'V,;'. 

power  spectral  density  function.  -  : 

For  the  detailed  discussion  and  the  estimation  of  the 
evolutionary  power  spectrum,  the  readers  are  referred  to  Refs.  18 
and  19. 

The  direct  generalization  of  the  above  discussion  to  n-dimen- 
sional  process  is  obvious.  Thus,  if  a  real  nonhomogeneous  process 
has  an  evolutionary  power  spectral  density  function,  the  process 
f^  (x)  can  be  simulated  by 

'  ■  '  ^  '  2  i';:-  ' 

fg  (x)  =  /2  S  [B  (x,cnj^)  S  A  03]  *  cos  (^  *  x+  f^)  (55) 

^  k=l 

where  B  (x,^)  is  the  n-dimensional  modulating  function  and  the 
remaining  notations  are  the  same  as  in  the  case  of  a  homogeneous 
process  given  by  Eq.  33.  It  can  be  shown  that  the  simulated  process 
(x)  possesses  the  target  evolutionary  power  spectrum  as 
N  -•  ‘®. 

Note  that  in  a  particular  case  when  B (x, m)  =  B  (x) ,  then 
f^ (x)  can  be  obtained  by  multiplying  a  homogeneous  process 
simulated  from  (m)  by  the  spatial  envelope  function  B (x) . 

A  more  detailed  study  with  numerical  examples  on  this  method 
of  simulation  for  nonhomogeneous  Gaussian  process  with  an  evolu- 
.  tionary  power  spectral  density  has  been  made  by  Yang  [6] . 
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In  the  following  section,  the  Gaussian  nonhomogeneous  " '  '  ; 

process  with  an  evolutionary  power  spectral  density  is  referred 
to  as  Gaussian  evolutionary  process  for  simplicity. 
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6.  Monte  Carlo  Solution  of  Structural  Dynamics 

■  v  Ths  preceding  method  of  digital  generation  of  sample  ^ 
functions  of  a  Gaussian  process  can  be  used  for  the  Monte  .. 
Carlo  solution  of  the  following  structural  problems.  It 
is  pointed  out  parenthetically  that  by  adding  a  constant 
value  m  to  the  sample  functions  f  (t,x)  described  in  the 
preceding  sections,  one  can  generate  sample  functions  of 
the  simulated  process  f  (t, x)  +  m  associated  with 
f,  (t,x)  +  m.  Note  that  the  mean  value  of  fa  (t,x)  +  m 
is  no  longer  zero  but  it  is  equal  to  m, 

(a)  The  method  can  be  used  in  the  response  analysis 
of  a  nonlinear  structure  under  random  loading  if  such  loading 
can  be  idealized  as  Gaussian  homogeneous  or  Gaussian  evo¬ 
lutionary  process  with  constant  mean  values.  In  particular, 
if  the  modes  of  the  corresponding  linear  structures 

are  known,  the  solution  (t,x)  is  in  approximation  ex¬ 
panded  into  a  finite  series. 


Yo  (t,x) 


K 

E 

k=l 


q^(t)  M^(x) 


(56) 


When  Eq.  56  is  substituted  into  the  governing  (nonlinear 
partial)  differential  equation  (s)  of  motion,  one  can  usually 
get  a  set  of  simultaneous  nonlinear  butordinary  differential 
equations  involving  the  generalized  forces  of  the 
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where  (t,x)  is  the  random  excitation  process  and  D 
indicates  an  appropriate  domain  of  integration.  Sample 
functions  then  be  digitally  generated  from 

Equation  57  with  fg  (t,x)  replaced  by  f(t;x); 

"  •I'  \  f  (t,  x)  dx  (58) 


It  goes  without  saying  that  f (t, x)  +  m  has  to  be  used 
in  place  of  f  (t,x)  if  the  excitation  process  has  a  non¬ 


zero  constant  mean  value  since  the  simple  superposition  of 


solutions  does  not  apply  in  this  case  because  of  nonlinearity.  ^ 
The  modes  M,  (x)  often  take  the  form  of  sinusoidal 


or  hyperbolic  functions  or  their  combinations.  Therefore, 


the  integration  in  Eq.  58  can  usually  be  carried  out  in 

!  ■  '  ■  ■  ■  '  '  '  '  ■  ■ '  ' 

closed  form  since  f  (t,x)  is  given  as  a  sum  of  cosine 


functions.  This  is  one  of  the  significant  advantages  of 
the  present  method  of  simulation  over  other  existing  methods. 


In  fact,  if  the  domain  of  integration  D  represents  a 


two  or  three  dimensional  space,  the  numerical  integration 


of  Eq.  58  will  usually  become  an  insurmountable  obstacle. 


Another  advantage  is  that  the  present  method  does  not 
require  the  nonlinearity  to  be  small  or  moderate,  a  condition 
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which  has  to  be  imposed  for  standard  linearization  or  per-r^^^^^;^ 
turbation  techniques.  ' 

Once  the  sample  functions  F  (t)  are  evaluated  from^  ^ 

Eq.  58,  then  the  sample  functions  q^  (t)  of  q^  (t)  can  be 
numerically  evaluated  from  the  (simultaneous)  nonlinear 

but  ordinary  differential  equations  mentioned  above  (replacing 
of  course  therein).  The  experience  shows 

that  this  phase  of  numerical  work  is  not  a  serious  problem. 
Finally,  the  sample  function  y(t,x)  of  the  solution 

Yo  {'t,x)  can  be  obtained  from  Eq.  56  with  replaced 

^  -■  '  -2  ■ 

t>Y  temporal  average  of  y  (t,x)  over  a  suffi¬ 

ciently  long  period  of  time  will  produce  the  mean  square 
response  in  the  Monte  Carlo  sense  if  the  processes  involved 
are  ergodic.  Otherwise,  the  ensemble  average  has  to  be 
considered. 


Reference  [20]  represents  a  typical  example  of  such 
analysis.  A  segment  of  a  sample  function  of  the  tip  deflection 
U(0,t)  of  a  vertical  pile  of  uniform  cross-section  in 
deep  water  (Fig.  6)  having  nonlinear  drag  effect  and  sub¬ 
jected  to  unidirectional  wind-induced  waves  is  shown  in 
Fig.  7,  where  a  segment  of  a  sample  function  of  the  response 
of  the  corresponding  linear  pile  (without  drag  term)  is 
also  shown  for  comparison.  In  this  study,  the  excitation 
is  due  to  waves  under  fully  developed  sea  conditions  with 
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mean  wind  velocity  V  for  which  the  Pierson  -  Moskowitz^ 
spectrum  (cu)  for  the  ocean-surface  elevation  has  been 
used;  ■  •  ' 


S"(m)  =  2SL  (59) 

where  a  =  8.10  x  10  ,  P  =  0.74  and  cDo  =  g/V  with 

g  s=  acceleration  due  to  gravity  and  V  =  mean  wind  velocity. 

The  application  of  this  type  of  Monte  Carlo  approach 
has  also  been  made  to  other  nonlinear  structural  response 
analysis  [3,4,21,22,23]. 


(b)  The  method  can  be  applied  to  the  failure  analysis 
of  a  structure  with  spatially  random  , variation  of  strength 
and  other  material  properties.  In  this  case,  sample  struc¬ 
tures  are  generated  by  digitally  generating  such  spatial 
variations  of  strength  and  other  material  properties. 

When  correlated  spatial  variations  are  observed  on  more 
than  one  material  property  (e.g.  Young's  modulus  and 
density),  usually  a  multidimensional  multivariate  process 
has  to  be  generated  with  the  aid  of  the  method  described 
in  Section  4. 

Applying  to  each  of  these  sample  structures  a  sample 
stress  history  of  a  random  stress  process,  the  fatigue 
life  of  a  sample  structure  can  be  computed  under  the 
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assumption  of  a  certain  fatigue  failure  mechanism; 
statistical  variation  of  the  fatigue  life  thus  computed^’^^^.^ 
establishes  its  empirical  distribution  under  the  random 
stress  process  in  the  Monte  Carlo  sense.  This  approach 
was  successfully  taken  in  Ref.  [24] .  A  similar  problem  in 
which  the  empirical  distribution  of  the  static  failure 

load  is  to  be  found  for  a  concrete  structure  with  spatial 
strength  variation  is  treated’ in  detail  in  Ref.  [25]. 

(c)  The  method  can  be  employed  effectively  when  the 
structural  system  to  be  analyized  is  complex  even  though 
it  involves  neither  nonlinearity  nor  random  variation  of 
material  properties.  The  mean  square  responses  (displace- 

.V  ■ 

ment,  shear  force  and  bending  moment)  of  a  large  floating 
plate  to  wind- induced  random  ocean-waves  are  computed  in 
Ref.  [26]  taking  the  temporal  averages  of  sample  response 
functions  as  in  Ref.  [20].  The  analysis  is  essentially 
numerical  since  sample  functions  of  the  wind-induced  ocean- 
surface  elevation  are  digitally  generated  and  the  corres¬ 
ponding  response  functions  are  numerically  obtained.  This 
was  done  because  the  ocean-structure  system  considered  was 
too  complex  to  solve  analytically.  Another  example  of  this 
kind  is  the  study  of  the  dynamic  interaction  between  moving 
vehicles  and  a  bridge  with  random  pavement  surface  rough- 
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ness  [27]  .  In  this  problem  the  random  pavement  surface  j 
roughness  is  digitally  simulated  for  numerical  response^  5^?^ 
analysis. 

In  some  problems  of  mechanics,  a  structure  is  considered  : 
complex  when  its  material  properties  are  spatially  random.  . 

The  wave  propagation  through  a  random  medium  is  one  of  these 
problems.  In  Ref.  [28] ,  the  stress  wave  propagation  through  ^ 

a  finite  cylinder  with  random  material  properties  is  treated 
under  the  condition  that  the  one  end  of  the  cylinder  is  acted 
upon  by  an  impact  load  and  the  other  end  is  free.  A  set 
of  one  hundred  samples  of  correlated  random  material  properties 
(Young's  modulus  and  density)  are  generated  thus  producing 
one  hundred  sample  cylinders.  The  finite  element  method  is 
applied  for  the  stress  analysis  to  compute  maximum  stress 
intensity  in  each  of  these  cylinders  due  to  the  impact.  An 
empirical  distribution  of  the  maximum  stress  intensity  is 
then  established  in  the  Monte  carlo  sense. 

(d)  Finally,  the  method  is  often  useful  when  the 
problem  is  to  determine  eigenvalues  (frequencies  and  buckling 
loads)  of  the  structure  with  random  material  properties.  As 
in  the  case  of  the  wave  propagation  problems,  sample  struc¬ 
tures  are  generated  and  the  statistical  distribution  of 
eigenvalues  of  these  structures  are  treated  as  the  empir  - 
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example  of  this  problem  is  given  in  Ref .  [29jv.s 
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Fig.  1  One-sided  spectral  density 


Pig.  2  Sample  functions  of  a  narrow-band  random  process  and 
its  envelope  process 

Fig.  3  Sample  functions  of  a  wide-band  random  process  and  its 
^envelope  process 

Fig.  4  Two-dimensional  spectral  density 

Pig^  5  Simulation  of  wind  velocity  at  different  points 


Pig.  6  Ocean  -  pile  system 


Fig.  7  A  section  of  sample  response  function  at  mean  wind 
velocity  23.6  ft/sec 
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Pig.  3  Sample  Functions  of  a  wide-band  random  process 
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Fig.  4  Two-dimensional  spectral  density 
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